## Modelo básico de hiv de Perelson

params <- c(
	c = 23,
	p = 700,  # ~
	lambda = 8,
	d = 0.01,
	delta = 1
)

a <- c(
	

	"dT1*x1",
	"drE*(x6/nu)*x2",
	"dT2*x3",
	"dT2*x4",
	"cV*x5",
	"(dE1+dmE*(x2/(x2+kd)))*x6",
	"dE2*x7",
	#####################
	"bT1*(x5/nu)*x1",
	"bT2*(x5/nu)*x3",
	#####################
	"gT*x1",
	"gT*x2",
	"gE*((x1+x2)/(x1+x2+kg))*x6",
	#####################
	"(aT*(x5/(x5+kV))+aA)*x3",
	"(aT*(x5/(x5+kV))+aA)*x4",
	"aE*(x5/(x5+kV))*x7",
	#####################
	"dV*x2",
	"(cT*(ks/(x5+ks)))*nu",
	"nu*cE+bE1*(x2/(x2+kb1))*x6",
	"bE2*(kb2/(x7+kb2))*x7"
)

V <- t(matrix(c(
	-1, 0, 0, 0, 0, 0, 0,
	 0,-1, 0, 0, 0, 0, 0,
	 0, 0,-1, 0, 0, 0, 0,
	 0, 0, 0,-1, 0, 0, 0,
	 0, 0, 0, 0,-1, 0, 0,
	 0, 0, 0, 0, 0,-1, 0,
	 0, 0, 0, 0, 0, 0,-1,
	#####################
	-1,+1, 0, 0,-1, 0, 0,
	 0, 0,-1,+1,-1, 0, 0,
	#####################
	-1, 0,+1, 0, 0, 0, 0,
	 0,-1, 0,+1, 0, 0, 0,
	 0, 0, 0, 0, 0,-1,+1,
	#####################
	params['nT']*+1,0,-1,0,0,0,0,
	0,params['nT']*+1,0,-1,0,0,0,
	0,0,0,0,0,params['nE']*+1,-1,
	#####################
	0,-1,0,0,params['nV']*+1,0,0,
	0, 0,+1, 0, 0, 0, 0,
	0, 0, 0, 0, 0,+1, 0,
	0, 0, 0, 0, 0, 0,+1
), ncol=7, byrow=TRUE))

fn_ode <- function(t, X, params) {
	with(as.list(c(X, params)), {
		# I use a double 'dd' prefix to avoid name collision
		ddT <- lambda - d*T - beta*V*T
		
		
		ddT1  <- - dT1*T1 - bT1*VI*T1 - gT*T1 + nT*((aT*VI)/(VI+kV)+aA)*T2
		ddT1_ <- bT1*VI*T1 - dV*T1_ - drE*E1*T1_ - gT*T1_ + nT*((aT*VI)/(VI+kV)+aA)*T2_
		ddT2  <- cT*(ks/(VI+ks)) + gT*T1 - dT2*T2 - bT2*VI*T2 - ((aT*VI)/(VI+kV)+aA)*T2
		ddT2_ <- gT*T1_ + bT2*VI*T2 - dT2*T2_ - ((aT*VI)/(VI+kV)+aA)*T2_
		ddVI  <- nV*dV*T1_ - cV*VI - (bT1*T1 + bT2*T2)*VI
		ddE1  <- cE + ((bE1*T1_)/(T1_+kb1))*E1 - ((dmE*T1_)/(T1_+kd))*E1 - dE1*E1 - gE*((T1+T1_)/(T1+T1_+kg))*E1 + nE*((aE*VI)/(VI+kV))*E2
		ddE2 <- gE*((T1+T1_)/(T1+T1_+kg))*E1 + ((bE2*kb2)/(E2+kb2))*E2 - dE2*E2 - ((aE*VI)/(VI+kV))*E2
		list(c(ddT1, ddT1_, ddT2, ddT2_, ddVI, ddE1, ddE2))
	})
}

X0 <- params['nu'] * c(x1=5, x2=1, x3=1400, x4=1, x5=10, x6=5, x7=1)

## algorithms

formatvector <- function(X) {
	# fit a given vector to a time vector of type: 0,1,2,3,...
	Y <- matrix(nrow=100+1, ncol=ncol(X))
	i <- 1
	for(t in seq(0,100,by=1)) {
		while(X[i,1] < t && i < nrow(X))
			i <- i+1
		Y[t+1,] <- c(t,X[i,-1])
	}
	colnames(Y) <- colnames(X)
	Y
}

myode <- function() {
	library(deSolve)
	X0_ <- X0
	names(X0_) <- c('T1','T1_','T2','T2_','VI','E1','E2')
	X <- ode(X0_, seq(0,100,by=0.01), fn_ode, params)
	colnames(X) <- c('time',paste('x',1:7,sep=''))  # convert back to old nomenclature
	formatvector(X)
}

mygillespie1 <- function(runid, X0, a, V, params, format) {
	library(GillespieSSA)
	X <- ssa(X0, a, V, params, 100, 'BTL')$data  # temp: time
	rownames(X) <- 1:nrow(X)
	colnames(X)[1] <- 'time'
	format(X)
}

mygillespies <- function(nruns) {
	library(parallel)
	cl <- makeCluster(detectCores())
	res <- parLapply(cl, 1:nruns, mygillespie1, X0=X0, a=a, V=V, params=params, format=formatvector)
	stopCluster(cl)
	## return standard devitions of these trajectories (ugly)
	#Y <- matrix(nrow=11, ncol=ncol(X))
	#for(t in 1:11) {
	#	Y[t,1] <- t-1
	#	for(x in 2:8) {
	#		v <- c()
	#		for(i in 1:nruns) {
	#			d <- res[[i]]
	#			v <- c(v,d[t,x])
	#		}
	#		Y[t,x] <- sd(v)
	#	}
	#}
	#colnames(Y) <- colnames(res[[1]])
	#Y
	res
}

## plot

library(gridExtra)
library(ggplot2)
library(reshape2)

theme_set(theme_gray(base_size=12) %+replace% theme(
	axis.text        = element_text(color='black'),
	panel.grid.minor = element_line(color='gray90'),
	panel.grid.major = element_line(color='gray60'),
	panel.background = element_rect(fill='white'),
	legend.text.align=0,
	legend.key = element_blank(),
	legend.position  = 'bottom'
))

X <- myode()
Y <- mygillespies(20)

#df <- melt(as.data.frame(X), id='time', variable.name='type', value.name='counts')

## add standard deviations to data.frame (ugly)
#s <- rep(0, nrow(df))
#for(i in 1:nrow(df))
#	s[i] <- Y[df[i,'time']+1, as.numeric(df[i,'type'])+1]
#df$s <- s

labels <- c(x1=expression(CD4[act]), x2=expression(CD4[act]^'\u2217'), x3=expression(CD4[rest]), x4=expression(CD4[rest]^'\u2217'), x5=expression(V[I]), x6=expression(CD8[eff]), x7=expression(CD8[mem]))

df <- melt(as.data.frame(X), id='time', variable.name='type', value.name='counts')
p <- ggplot()
p <- p + geom_line(aes(x=time, y=counts, color=type, size=type, linetype=type), df)

for(y in Y) {
	df <- melt(as.data.frame(y), id='time', variable.name='type', value.name='counts')
	p <- p + geom_line(aes(x=time, y=counts, color=type, size=type, linetype=type), df, alpha=0.2)
}

p <- p + scale_y_log10()
p <- p + scale_colour_manual(name="Type", values=c(x1='blue',x2='blue',x3='blue',x4='blue', x5='red', x6='green',x7='green'), labels=labels)
p <- p + scale_fill_manual(name="Type", values=c(x1='blue',x2='blue',x3='blue',x4='blue', x5='red', x6='green',x7='green'), labels=labels)
p <- p + scale_size_manual(name="Type", values=c(x1=1,x2=1,x3=.5,x4=.5,x5=.5,x6=1,x7=.5), labels=labels)
p <- p + scale_linetype_manual(name="Type", values=c(x1=1,x2=2,x3=1,x4=2,x5=1,x6=1,x7=1), labels=labels)
print(p)
ggsave('plot.pdf', p, width=10, height=5)
